Hi,
I want to have a colorbar aside the curve, but I can't. (My Maple version is 2018).

restart:
with(plots):

# --- constants ---
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# --- parameter ranges ---
deltac_min := 0.01: deltac_max := 0.99:
Tch_min := 0.1:     Tch_max := 1.1:

# --- numerical safety ---
epsDen := 1e-12:
Mcap := 3.0: # cap to avoid infinity

# --- denominator ---
Den := (deltac, Tch) ->
       -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
       + 2*deltac*qc + 2*Tch + 2*deltac:

# --- safe functions ---
M1safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

M0safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2*alpha/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

# --- density plots ---
p1 := densityplot(
       M1safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M1(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

p2 := densityplot(
       M0safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M0(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

# --- Contour Den = 0 (to mark physical boundary) ---
cont := contourplot(
          Den(deltac,Tch),
          deltac = deltac_min .. deltac_max,
          Tch   = Tch_min .. Tch_max,
          contours = [0],
          color = black,
          thickness = 2
        ):

# --- Overlay contour ---
p1_final := display(p1, cont):
p2_final := display(p2, cont):

# --- Display side by side ---
display(array([p1_final, p2_final]), scaling = constrained);

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

restart;
with(plots):

# === Constants ===
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# === Denominator ===
Den := (deltac, Tch) ->
     -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
     + 2*deltac*qc + 2*Tch + 2*deltac:

# === M1 and M0 ===
M1 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2/sqrt(Den(deltac,Tch)), undefined):
M0 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2*alpha/sqrt(Den(deltac,Tch)), undefined):

# === Domain ===
deltac_min := 0.05:  deltac_max := 0.95:
Tch_min := 0.01:     Tch_max := 1.10:   # Physically realistic range

# === 3D Plots ===
p1 := plot3d(M1(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=green,
             labels=["δ_c", "T_ch", "M1"],
             labelfont=[TIMES,12],
             title="M1(δ_c, T_ch) for qc=0.4, qh=0.64"):

p2 := plot3d(M0(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=blue,
             labels=["δ_c", "T_ch", "M0"],
             labelfont=[TIMES,12],
             title="M0(δ_c, T_ch) for α=0.8, qc=0.4, qh=0.64"):

# === Display vertically (M1 above M0) ===
display(Array(1..2, [p1, p2]), insequence=true, scaling=constrained);
display(Array([[p1], [p2]]));

Good day.

Using the iterative map routine, the location of the bifurcation points of the logistic function can be determined using the plot (see attached).

I was wondering .. is it possible to estimate and output the locations of these points and the range of the function? I would like to explore the bifurcation behavior for the modified function so, this would be a great help. 

In this case, the location of the points are: (1.00,0.00), (3.00, 0.66), (3.45, 0.44), (3.45, 0.85) and the range is [0,4]. 

Thanks for reading!

MaplePrimes_Oct_16.mw

Good day.

I am looking into the behaviour of a function, V, that depends on several parameter values; these values are fixed in the attached example. However, I have encountered an issue that is puzzling me and I was hoping that someone may be able to shine a light on this for me.

Basically,  I would like to understand how the solution, V, behaves as the value of exponent, beta, approaches infinity.

Straightforward analysis suggests that the value of V tends to -10 as beta grows infinitely large (s -> C, and beta ->infinity and so, V -> -10 ... so far, so good).

However, when the function is plotted, the solution seems to converge to a value, 6.5, as beta tends to a very large number (10^17).

Now .. here's the mystery .. there appears to be a critical value of around 7.854 x 10^17; here, the limit seems to switch from V=6.5 to -10. Does this phenomenon correspond to a discontinuity or is it related to the computational process? Are there any built-in routines in Maple to check for such potential conditions?

Thanks for reading!

MaplePrimes_Sep_19.mw

dear sir here i am giving python code exicuting 3d plots but i cant plot animations like

restart;
Nc := .3; M := 1.5; QG := 1.5; Thetaa := .2; n1 := 1; n2 := 0; lambda1 := .1; lambda2 := .1; lambda3 := .1;

p := 2; a := .5; alpha1 := (1/2)*Pi;

p1 := 0.1e-1; p2 := 0.1e-1;
rf := 997.1; kf := .613; cpf := 4179; betaf := 21*10^(-5);

betas1 := .85*10^(-5); rs1 := 3970; ks1 := 40; cps1 := 765;
betas2 := 1.67*10^(-5); rs2 := 8933; ks2 := 401; cps2 := 385;

z1 := 1/((1-p1)^2.5*(1-p2)^2.5);
knf := kf*(ks1+2*kf-2*p1*(kf-ks1))/(ks1+2*kf+p1*(kf-ks1)); khnf := knf*(ks2+2*knf-2*p2*(knf-ks2))/(ks2+2*knf+p2*(knf-ks2));
z2 := 1-p1-p2+p1*rs1/rf+p2*rs2/rf;
z3 := 1-p1-p2+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf);
z4 := khnf/kf;
z5 := 1-p1-p2+p1*rs1*betas1/(rf*betaf)+p2*rs2*betas2/(rf*betaf);
OdeSys := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Cond := {Theta(0, tau) = 1+a*sin(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSys1 := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Cond1 := {Theta(0, tau) = 1+a*sin(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSysa := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Conda := {Theta(0, tau) = 1+a*cos(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSysa1 := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Conda1 := {Theta(0, tau) = 1+a*cos(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};

colour := [cyan, black];
colour1 := [red, blue];
NrVals := [2.5, 3.5];
for j to numelems(NrVals) do Ans := pdsolve((eval([OdeSys, Cond], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plots[j] := Ans:-plot(Theta(X, tau), tau = .8, linestyle = "solid", labels = ["Y", Theta(Y, tau)], color = colour[j], 'axes' = 'boxed'); eta[j] := Ans:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n1*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "solid", 'axes' = 'boxed', labels = [" τ ", " η "], color = colour[j]); Ans1 := pdsolve((eval([OdeSys1, Cond1], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plots1[j] := Ans1:-plot(Theta(X, tau), tau = .8, linestyle = "dash", labels = ["Y", Theta(Y, tau)], color = colour[j], 'axes' = 'boxed'); eta1[j] := Ans1:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n2*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "dash", 'axes' = 'boxed', labels = [" τ ", " η "], color = colour[j]); Ansa := pdsolve((eval([OdeSysa, Conda], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plotsa[j] := Ansa:-plot(Theta(X, tau), tau = .8, linestyle = "solid", labels = ["Y", Theta(Y, tau)], color = colour1[j], 'axes' = 'boxed'); etaa[j] := Ansa:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n1*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "solid", 'axes' = 'boxed', labels = [" τ ", " η "], color = colour1[j]); Ansa1 := pdsolve((eval([OdeSysa1, Conda1], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plotsa1[j] := Ansa1:-plot(Theta(X, tau), tau = .8, linestyle = "dash", labels = ["Y", Theta(Y, tau)], color = colour1[j], 'axes' = 'boxed'); etaa1[j] := Ansa1:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n2*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "dash", 'axes' = 'boxed', labels = [" τ ", " η "], color = colour1[j]) end do;
plotA := plots:-display([seq(Plots[j], j = 1 .. 2)]);
plotB := plots:-display([seq(Plots1[j], j = 1 .. 2)]);
plotAA := plots:-display([seq(Plotsa[j], j = 1 .. 2)]);
plotBB := plots:-display([seq(Plotsa1[j], j = 1 .. 2)]);
plotC := plots:-display([seq(eta[j], j = 1 .. 2)]);
plotD := plots:-display([seq(eta1[j], j = 1 .. 2)]);
plotCC := plots:-display([seq(etaa[j], j = 1 .. 2)]);

plotDD := plots:-display([seq(etaa1[j], j = 1 .. 2)]);
plots:-display([plotA, plotB, plotAA, plotBB], size = [700, 700]);

plots:-display([plotC, plotD, plotCC, plotDD], size = [700, 700]);

how to take animation of the plots

i have seen some plots in maple also for that reason i have posted this question here

paper2_new_efficiency_plots_2025.mw

Hi,
I want to solve equations (1) and (2) to obtain the values of certain quantities, for example x and f. In fact, both equations must be satisfied simultaneously for \tau and f. Since there are other parameters in these equations (for example \tau, M, etc.), and I do not know the optimal values of these parameters, I would like to solve equations (1) and (2) over given ranges for tau, M, and so on. I know that I should use loop commands, but I keep encountering errors. Please guide me.
fsolve.mw

Hi 
How can I plot f3=0 for different values of Tch and qc?

1.mw

How to get amultiple line plots in one graph for differnt values of E1 := .1,0.2,0.3,and 0.4 with differnt color line red,blue,black and green and how to get numerical values for the all E1 value in one matrix form 

at E1=0.1 value of diff(g(x), x, x), diff(f(x), x, x),diff(f(x), x) and diff(g(x),  x).

AGM_method_single_line_plot.mw

How to solve two boundary problems in one graph not getting graphs shown in pdf
symetry_paper_work.mw
symmetry_graphs_pdf.pdf

How to get table values 

why this error is getting. but in the published paper for the same parameter it is converging. what is the mistake in this worksheet please rectify it in sachin_base_paper.mw

Someone should please help me compute the left and right eigenvectors of the system below. The purpose is to compute values for 'a' and 'b' in the bifurcation formula.

Thank you

Download CBD2.mw

It is not difficult to manually check the validity of the identity  sqrt(x+2*sqrt(x-1)) = sqrt(x-1)+1 , which is true at least for all  x>=1 . I don't know of any Maple command that directly makes this simplification:

expr:=sqrt(x+2*sqrt(x-1));
simplify(expr);
simplify(expr) assuming x>=1;

It was possible to simplify it only by making a variable substitution  x=y+1 (and then a reverse substitution  y=x-1 ) :

expr1:=subs(x=y+1,expr);
simplify(expr1);
subs(y=x-1, %);

By the way,  the CAS Mathematica also cannot cope with simplifying the expression  expr .

Good day.

I am constructing a 4-set Venn Diagram and I would like to know if it is possible to fix the number of decimal places in the solution.

The attached worksheet is given as an example; the default number of decimal places seems to be 2. I would like this to be either 0 or 1 (for both absolute and relative values). 

Does anyone know how to do this? 

Thanks for reading!

MaplePrimes_Venn_Diagram.mw

How to solve these pde equations in maple to get the similar type graphs.

Ode equations we can solve directly but these equations are pde .

in the article they have solved the governing equations by series solution? 

can we solve these equations in maple also by series solution or any other method is there to solve these equations

Someone please help me with the computation of the right and left eigenvectors. my system of equation is attached below

Download Cotton_DFE_and_Jacobian.mw

Please, I am encountering error trying to run these codes for sensitivity analysis using the formula for sensitivity analysis

Download Computed_Sensitivity_Analys_for_CBD.mw